/********************************
 * Robinson Matrix Algorithm
 *
 * Copyright (C) 2009 CNRS
 * Author : Florent AVELLANEDA, Eugen CEPOI
 * Algorithm design : Victor CHEPOI, Morgan SESTON ( http://www.lif-sud.univ-mrs.fr/%7Echepoi/robinson_chepoi_seston.pdf )
 * 
 *
 * All rights reserved.
 *
 *   This file is part of Robinson Matrix Algorithm.
 *
 *   Robinson Matrix Algorithm is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   Robinson Matrix Algorithm is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with Robinson Matrix Algorithm.  If not, see <http://www.gnu.org/licenses/>.
 *
 ***********************************************************************************/
#include "Graphs.h"

#include <assert.h>

#include "Matrix.h"
#include "PartialOrder.h"
#include "Util.h"
#include "AlgoGraph.h"
#include "Var.h"
#include "SAT.h"

using namespace Util;

namespace Graphs
{
    bool partitionAndSort(Matrix *m, double epsilon, unsigned int leftBoundHole1, unsigned int rightBoundHole1,
            unsigned int leftBoundHole2, unsigned int rightBoundHole2, std::set<unsigned int> &Xij,
            std::vector<std::set<unsigned int> > &cellsOfOL, std::set<unsigned int> &X_m,
            std::set<unsigned int> &X_p, std::vector<PartialOrder*> &Cij_m,
            std::vector<PartialOrder*> &Cij_p, std::vector<unsigned int> &cellsOrder_m,
            std::vector<unsigned int> &cellsOrder_p)
    {
        // Define struct {{{
        std::vector<double> d;

        std::vector<std::set<unsigned int> > S;

        std::vector<std::set<unsigned int> > L;
        std::vector<unsigned int> blockOfL;

        std::vector<std::set<unsigned int> > OL;
        std::vector<unsigned int> xCellsOfOl;

        std::vector<std::set<unsigned int> > psy;
        std::vector<unsigned int > clusterOfPsy;

        std::vector<std::set<unsigned int> > Gij;
        std::vector<std::set<unsigned int> > G1;
        std::vector<std::set<unsigned int> > G2;
        std::vector<std::set<unsigned int> > G3;
        // }}} End define struct

        std::vector<unsigned int> X;
        {
            for (std::set<unsigned int>::iterator it = Xij.begin(); it != Xij.end(); ++it)
                X.push_back(*it);
        }

        calculateDx(d, X, m);
        graphS(S, X, d, epsilon, m);
        graphL(L, X, d, epsilon, m);
        partitionL(blockOfL, L);
        graphOL(OL, X, blockOfL, L, d, epsilon, m);
        cellsOL(L, OL, cellsOfOL, xCellsOfOl);
        graphPsy(psy, cellsOfOL, xCellsOfOl, blockOfL, S, L, X, d, epsilon, m);

        //Testing if psy is biparti and getting the clusters
        std::vector<std::set<unsigned int> > clusters;
        if ( !AlgoGraph::bipartiGraphClusters(psy, clusters, clusterOfPsy) )
            return false;

        std::vector< SAT::Var > vVar;
        if( !sat(vVar, Gij, G1, G2, G3, cellsOfOL, clusterOfPsy, clusters, X, L, d, epsilon, m, leftBoundHole1, rightBoundHole1, leftBoundHole2, rightBoundHole2) )
            return false;

        partition(X_m, X_p, Gij, Cij_m, Cij_p, leftBoundHole1, rightBoundHole1, leftBoundHole2,
                rightBoundHole2, vVar, cellsOrder_m, cellsOrder_p, cellsOfOL, X, L, d,
                epsilon, m);

        return true;
    }

    void graphS(std::vector<std::set<unsigned int> > &S, std::vector<unsigned int> &X,
            std::vector<double> &d, double epsilon, Matrix *m)
    {
        S.resize(X.size());
        for (unsigned int i = 0; i < X.size(); i++)
            for (unsigned int j = i+1; j < X.size(); j++)
            {
                if( m->get(X[i], X[j]) > (max( d[i], d[j] ) + 3.0*epsilon) )
                {
                    S[i].insert(j);
                    S[j].insert(i);
                }
            }
    }

    void graphL(std::vector<std::set<unsigned int> > &L, std::vector<unsigned int> &X,
            std::vector<double> &d, double epsilon, Matrix *m)
    {
        L.resize(X.size());
        for (unsigned int i = 0; i < X.size(); i++)
            for (unsigned int j = i+1; j < X.size(); j++)
            {
                if( m->get(X[i], X[j]) < (max( d[i], d[j] ) - 3.0*epsilon) )
                {
                    L[i].insert(j);
                    L[j].insert(i);
                }
            }
    }

    void partitionL(std::vector<unsigned int> &blockOfL, std::vector<std::set<unsigned int> > &L)
    {
        blockOfL.resize( L.size() );
        std::set<unsigned int> S;
        for (unsigned int i = 0; i < L.size(); i++)
            S.insert(i);
        AlgoGraph::connected<unsigned int>( L, S, blockOfL );
    }

    void graphOL(std::vector<std::set<unsigned int> > &OL, std::vector<unsigned int> &X,
            std::vector<unsigned int> &blockOfL, std::vector<std::set<unsigned int> > &L,
            std::vector<double> &d, double epsilon, Matrix *m )
    {
        OL.resize(X.size());
        for (unsigned int i = 0; i < X.size(); i++)
            for (unsigned int j = 0; j < X.size(); j++)
            {
                if(i==j)
                    continue;

                if( blockOfL[i] == blockOfL[j] )
                {
                    if( aarrow( i, j, epsilon, d, X, L, m ) )
                    {
                        OL[i].insert(j);
                        continue;
                    }
                }

                if( m->get(X[i], X[j]) < (max( d[i], d[j] ) - 5.0*epsilon) )
                    OL[i].insert(j);
            }
    }

    void cellsOL(std::vector<std::set<unsigned int> > &L, std::vector<std::set<unsigned int> > &OL,
            std::vector<std::set<unsigned int> > &cellsOfOL, std::vector<unsigned int> &xCellsOfOl )
    {
        std::set<unsigned int> Som;
        for (unsigned int i = 0; i < L.size(); i++)
            Som.insert(i);

        AlgoGraph::stronglyConnected<unsigned int, std::vector<std::set<unsigned int> > >(OL, Som, xCellsOfOl, cellsOfOL);
    }

    void graphPsy(std::vector<std::set<unsigned int> > &psy, std::vector<std::set<unsigned int> > &cellsOfOL,
            std::vector<unsigned int> xCellsOfOl, std::vector<unsigned int> &/*blockOfL*/,
            std::vector<std::set<unsigned int> > &S, std::vector<std::set<unsigned int> > &L,
            std::vector<unsigned int> &X, std::vector<double> &d, double epsilon, Matrix *m )
    {
        psy.resize( cellsOfOL.size() );

        // S1
        for (unsigned int i = 0; i < S.size(); i++)
            for (std::set<unsigned int>::iterator it = S[i].begin(); it != S[i].end(); ++it)
                if( xCellsOfOl[i] != xCellsOfOl[*it] )
                {
                    psy[ xCellsOfOl[i]   ].insert( xCellsOfOl[*it] );
                    psy[ xCellsOfOl[*it] ].insert( xCellsOfOl[i]   );
                }

        // S2 : TODO : Optimisation possible
        for (unsigned int i = 0; i < cellsOfOL.size(); i++)
        {
            for (std::set<unsigned int>::iterator it = cellsOfOL[i].begin(); it != cellsOfOL[i].end(); ++it)
            {
                std::set<unsigned int>::iterator it2 = it;
                ++it2;
                for ( ; it2 != cellsOfOL[i].end(); ++it2)
                {



                    for (unsigned int ip = i+1; ip < cellsOfOL.size(); ip++)
                    {
                        for (std::set<unsigned int>::iterator itp = cellsOfOL[ip].begin(); itp != cellsOfOL[ip].end(); ++itp)
                        {
                            std::set<unsigned int>::iterator it2p = itp;
                            ++it2p;
                            for ( ; it2p != cellsOfOL[ip].end(); ++it2p)
                            {

                                if( (aarrow( *it, *itp, epsilon, d, X, L, m) && aarrow( *it2, *it2p, epsilon, d, X, L, m )) ||
                                        (aarrow( *it, *it2p, epsilon, d, X, L, m) && aarrow( *it2, *itp, epsilon, d, X, L, m )) ||
                                        (aarrow( *itp, *it, epsilon, d, X, L, m) && aarrow( *it2p, *it2, epsilon, d, X, L, m )) ||
                                        (aarrow( *itp, *it2, epsilon, d, X, L, m) && aarrow( *it2p, *it, epsilon, d, X, L, m )) )
                                {
                                    psy[i].insert(ip);
                                    psy[ip].insert(i);
                                }

                            }
                        }
                    }
                }
            }
        }
    }

    bool isG3(unsigned int i, unsigned int j, std::vector<std::set<unsigned int> > &cellsOfOL,
            double epsilon, std::vector<unsigned int> &X,
            Matrix *m, const unsigned int boundA, const unsigned int boundB)
    {
        for (std::set<unsigned int>::iterator it = cellsOfOL[i].begin(); it != cellsOfOL[i].end(); ++it)
            for (std::set<unsigned int>::iterator it2 = cellsOfOL[j].begin(); it2 != cellsOfOL[j].end(); ++it2)
            {
                if( m->get(X[*it], boundA) < (m->get(X[*it2], boundA) - 16.0*epsilon) ||
                        m->get(X[*it], boundB) < (m->get(X[*it2], boundB) - 16.0*epsilon) )
                    return true;
            }

        return false;
    }


    bool sat( std::vector< SAT::Var > &vVar, std::vector<std::set<unsigned int> > &Gij, std::vector<std::set<unsigned int> > &G1, std::vector<std::set<unsigned int> > &G2,
            std::vector<std::set<unsigned int> > &G3, std::vector<std::set<unsigned int> > &cellsOfOL,
            std::vector<unsigned int> &clusterOfPsy, std::vector<std::set<unsigned int> > &clusters,
            std::vector<unsigned int> &X, std::vector<std::set<unsigned int> > &L,
            std::vector<double> &d, double epsilon, Matrix *m, unsigned int leftBoundHole1,
            unsigned int rightBoundHole1, unsigned int leftBoundHole2, unsigned int rightBoundHole2)
    {
        // Generate graph G {{{
        G1.resize( cellsOfOL.size() );
        G2.resize( cellsOfOL.size() );
        G3.resize( cellsOfOL.size() );
        Gij.resize( cellsOfOL.size() );
        for (unsigned int i = 0; i < cellsOfOL.size(); i++)
            for (unsigned int j = 0; j < cellsOfOL.size(); j++)
            {
                if(i==j)
                    continue;
                if( twinClusters( clusterOfPsy, i, j ) )
                {
                    G1[i].insert(j);
                    Gij[i].insert(j);
                }
                else if( isG2(i, j, cellsOfOL, d, epsilon) )
                {
                    Gij[i].insert(j);
                    G2[i].insert(j);
                }
                else if( isG3(i, j, cellsOfOL, d, epsilon, X, L, m) )
                {
                    Gij[i].insert(j);
                    G3[i].insert(j);
                }
            }
        // }}} End generate graph G

        // Test G3 Clycle {{{
        std::vector<unsigned int> maximalChain;
        if( !isAcyclique(G3) )
            return false;
        // }}} End test G3 clycle

        // Generate 2-SAT problem {{{
        // vVar[ indice_de_cellsOfOL ].getVal() = X- (false) ou X+ (true)
        vVar.resize( cellsOfOL.size() );
        SAT::SAT sat;

        // F1
        for (unsigned int i = 0; i < clusters.size(); i++)
        {
            for (std::set<unsigned int>::iterator it = clusters[i].begin(); it != clusters[i].end(); ++it)
            {
                std::set<unsigned int>::iterator it2 = it;
                for (++it2; it2 != clusters[i].end(); ++it2)
                {
                    sat.addClause( !vVar[ *it ] + vVar[ *it2 ] );
                    sat.addClause( vVar[ *it ] + !vVar[ *it2 ] );
                }
            }
        }

        // F2
        for (unsigned int i = 0; i < clusters.size(); i=i+2)
            for (std::set<unsigned int>::iterator it = clusters[i].begin(); it != clusters[i].end(); ++it)
                for (std::set<unsigned int>::iterator it2 = clusters[i+1].begin(); it2 != clusters[i+1].end(); ++it2)
                {
                    sat.addClause( vVar[ *it ] + vVar[ *it2 ] );
                    sat.addClause( !vVar[ *it ] + !vVar[ *it2 ] );
                }
        std::set< std::pair<unsigned int, unsigned int> > omega;

        createOmega( omega, G1, G2, G3 );

        // F3 & F4
        for (std::set< std::pair<unsigned int, unsigned int> >::iterator it = omega.begin(); it != omega.end(); ++it)
        {
            sat.addClause( vVar[ it->first ] + vVar[ it->second ] );
            sat.addClause( !vVar[ it->first ] + !vVar[ it->second ] );
        }
        // }}} End Generate 2-SAT problem

        if ( !sat.solve2Sat() )
            return false;

        // find new G3 arc with bound
        unsigned int boundA, boundB;
        for (unsigned int i = 0; i < cellsOfOL.size(); i++)
            for (unsigned int j = 0; j < cellsOfOL.size(); j++)
            {
                if(i==j)
                    continue;

                if( Gij[i].count(j) == 0)
                {
                    if( vVar[i].getVal() == vVar[j].getVal() )
                    {
                        if( vVar[ i ].getVal() < 1 )
                        {
                            boundA = leftBoundHole1;
                            boundB = rightBoundHole1;
                        }else
                        {
                            boundA = leftBoundHole2;
                            boundB = rightBoundHole2;
                        }

                        if( isG3( i, j, cellsOfOL, epsilon, X, m, boundA, boundB ) )
                            Gij[i].insert(j);
                    }
                }
            }




        return true;
    }

    bool isG2(unsigned int i, unsigned int j, std::vector<std::set<unsigned int> > &cellsOfOL, std::vector<double> &d, double epsilon)
    {
        for (std::set<unsigned int>::iterator it = cellsOfOL[i].begin(); it != cellsOfOL[i].end(); ++it)
            for (std::set<unsigned int>::iterator it2 = cellsOfOL[j].begin(); it2 != cellsOfOL[j].end(); ++it2)
                if( d[ *it ] < (d[ *it2 ] - 4.0*epsilon) )
                    return true;

        return false;
    }

    bool isG3(unsigned int i, unsigned int j, std::vector<std::set<unsigned int> > &cellsOfOL,
            std::vector<double> &/*d*/, double epsilon, std::vector<unsigned int> &X,
            std::vector<std::set<unsigned int> > &L, Matrix *m)
    {
        for (std::set<unsigned int>::iterator it = cellsOfOL[i].begin(); it != cellsOfOL[i].end(); ++it)
            for (std::set<unsigned int>::iterator it2 = cellsOfOL[j].begin(); it2 != cellsOfOL[j].end(); ++it2)
                for (unsigned int z = 0; z < X.size(); z++)
                    if(     (z != *it) && (z != *it2) &&
                            (m->get(X[*it], X[z]) < (m->get(X[*it2], X[z]) - 16.0*epsilon)) &&
                            (L[*it].count( z ) == 0) &&
                            (L[*it2].count( z ) == 0)  )
                        return true;
        return false;
    }


    void createOmega( std::set< std::pair<unsigned int, unsigned int> > &omega,
            std::vector<std::set<unsigned int> > &G1, std::vector<std::set<unsigned int> > &G2,
            std::vector<std::set<unsigned int> > &G3 )
    {
        std::vector<std::set<unsigned int> > NG3;
        NG3.resize( G3.size() );

        for (unsigned int i = 0; i < G3.size(); i++)
            for (std::set<unsigned int>::iterator it = G3[i].begin(); it != G3[i].end(); ++it)
                NG3[ *it ].insert( i );

        for (unsigned int i = 0; i < G1.size(); i++)
            for (std::set<unsigned int>::iterator it = G1[i].begin(); it != G1[i].end(); ++it)
            {
                if(i == *it)
                    continue;
                omega.insert( std::pair<unsigned int, unsigned int>( i, *it) );
            }

        for (unsigned int i = 0; i < G2.size(); i++)
        {
            // Noeud a ne pas prendre en compte pour la recherche
            std::set<unsigned int> N1;
            for (std::set<unsigned int>::iterator it = G3[i].begin(); it != G3[i].end(); ++it)
                N1.insert( *it );

            bool omega2 = false;
            for (std::set<unsigned int>::iterator it = G2[i].begin(); it != G2[i].end(); ++it)
            {
                if(*it==i)
                    continue;

                // Noeud a ne pas prendre en compte pour la recherche
                std::set<unsigned int> N2 = N1;
                for (std::set<unsigned int>::iterator it2 = NG3[*it].begin(); it2 != NG3[*it].end(); ++it2)
                    N2.insert( *it2 );

                // Chercher omega3
                for (std::set<unsigned int>::iterator it2 = G2[*it].begin(); it2 != G2[*it].end(); ++it2)
                {
                    if(*it2 == *it)
                        continue;

                    // Noeud a ne pas prendre en compte pour la recherche
                    std::set<unsigned int> N3 = N2;
                    for (std::set<unsigned int>::iterator it3 = G3[*it].begin(); it3 != G3[*it].end(); ++it3)
                        N3.insert( *it3 );
                    for (std::set<unsigned int>::iterator it3 = NG3[*it2].begin(); it3 != NG3[*it2].end(); ++it3)
                        N3.insert( *it3 );

                    if( AlgoGraph::existsPath( G3, *it2, i, N3 ) )
                    {
                        omega2 = true;
                        omega.insert( std::pair<unsigned int, unsigned int>( *it, *it2 ) ); // OMEGA 3
                    }
                }

                // Chercher omega1
                if( AlgoGraph::existsPath( G3, *it, i, N2 ) )
                    omega.insert( std::pair<unsigned int, unsigned int>( i, *it ) ); // OMEGA 1
                else if( omega2 )
                    omega.insert( std::pair<unsigned int, unsigned int>( i, *it ) ); // OMEGA 2
            }
        }
    }


    void partition(std::set<unsigned int> &X_m, std::set<unsigned int> &X_p, std::vector<std::set<unsigned int> > &Gij,
            std::vector<PartialOrder* > &Cij_m, std::vector<PartialOrder* > &Cij_p, unsigned int leftBoundHole1,
            unsigned int rightBoundHole1, unsigned int leftBoundHole2, unsigned int rightBoundHole2,
            std::vector< SAT::Var > &vVar, std::vector<unsigned int> &cellsOrder_m,
            std::vector<unsigned int> &cellsOrder_p, std::vector<std::set<unsigned int> > &cellsOfOL,
            std::vector<unsigned int> &X, std::vector<std::set<unsigned int> > &L,
            std::vector<double> &d, double epsilon, Matrix *m )
    {
        std::vector<std::set<unsigned int> >Gij_m, Gij_p;
        Gij_m.resize(Gij.size());
        Gij_p.resize(Gij.size());

        // Partitioning Xij into Xij- and Xij+
        for ( unsigned int i = 0; i < vVar.size(); ++i )
            if ( vVar[i].getVal() < 1 )
            {
                X_m.insert(i);
                Gij_m[i] = Gij[i];
            }
            else
            {
                X_p.insert(i);
                Gij_p[i] = Gij[i];
            }
        // end Partitioning

        // Finding a topological order on the cells of Xij- and Xij+
        for ( std::set<unsigned int>::iterator i = X_m.begin(); i != X_m.end(); ++i )
            for ( std::set<unsigned int>::iterator j = X_p.begin(); j != X_p.end(); ++j )
                Gij_m[*i].erase(*j);


        for ( std::set<unsigned int>::iterator i = X_p.begin(); i != X_p.end(); ++i )
            for ( std::set<unsigned int>::iterator j = X_m.begin(); j != X_m.end(); ++j )
                Gij_p[*i].erase(*j);

        AlgoGraph::topologicalSorting(Gij_m, X_m, cellsOrder_m);
        AlgoGraph::topologicalSorting(Gij_p, X_p, cellsOrder_p);
        // End topological order


        /* Construction of the partial orders <=Cij- et <=Cij+ for each cell of Xij- and Xij+ */
        for ( unsigned int i = 0; i < cellsOrder_m.size(); ++i )
        {
            Matrix m_m(m, Util::transformer(cellsOfOL[cellsOrder_m[i]], X) );
            m_m.setElementFromMainMatrix(leftBoundHole1);
            m_m.setElementFromMainMatrix(rightBoundHole1);

            Cij_m.push_back(new PartialOrder(m_m));
            Cij_m.back()->setEpsilon( epsilon );
            bool atLeastOneRuleAdded = false;
            unsigned int pos1, pos2;
            pos1 = 0;

            for ( std::set<unsigned int>::iterator it = cellsOfOL[cellsOrder_m[i]].begin(); it != cellsOfOL[cellsOrder_m[i]].end(); ++it )
            {
                pos2 = pos1+1;
                std::set<unsigned int>::iterator it2 = it;
                ++it2;
                for ( ; it2 != cellsOfOL[cellsOrder_m[i]].end(); ++it2, ++pos2 )
                    if ( aarrow((*it), (*it2), epsilon, d, X, L, m) )
                    {
                        atLeastOneRuleAdded = true;
                        Cij_m.back()->setOrientation(pos1, pos2);
                    }
                ++pos1;
            }

            Cij_m.back()->setOrientation(m_m.getSize()-2, m_m.getSize()-1);
            Cij_m.back()->update( false );
            Cij_m.back()->setConstructBeetweness(!atLeastOneRuleAdded);
        }

        for ( unsigned int i = 0; i < cellsOrder_p.size(); ++i )
        {
            Matrix m_p(m, Util::transformer(cellsOfOL[cellsOrder_p[i]], X));
            m_p.setElementFromMainMatrix(leftBoundHole2);
            m_p.setElementFromMainMatrix(rightBoundHole2);

            Cij_p.push_back(new PartialOrder(m_p));
            Cij_p.back()->setEpsilon( epsilon );

            bool atLeastOneRuleAdded = false;
            unsigned int pos1, pos2;
            pos1 = 0;
            for ( std::set<unsigned int>::iterator it = cellsOfOL[cellsOrder_p[i]].begin(); it != cellsOfOL[cellsOrder_p[i]].end(); ++it )
            {
                pos2 = pos1+1;
                std::set<unsigned int>::iterator it2 = it;
                ++it2;
                bool atLeastOneRuleAdded = false;

                for ( ; it2 != cellsOfOL[cellsOrder_p[i]].end(); ++pos2, ++it2 )
                    if ( aarrow((*it), (*it2), epsilon, d, X, L, m) )
                    {
                        atLeastOneRuleAdded = true;
                        Cij_p.back()->setOrientation(pos1, pos2);
                    }
                ++pos1;
            }

            Cij_p.back()->setOrientation(m_p.getSize()-2, m_p.getSize()-1);
            Cij_p.back()->update( false );
            Cij_p.back()->setConstructBeetweness(!atLeastOneRuleAdded);
        }
    }

    bool aarrow(unsigned int i, unsigned int j, double epsilon, std::vector<double> &d,
            std::vector<unsigned int> &X, std::vector<std::set<unsigned int> > &L, Matrix *m)
    {
        if( d[i] < (d[j] - 4.0*epsilon) )
            return true;

        if( d[i] >= (d[j] - 4.0*epsilon) )
        {
            for (unsigned int k = 0; k < X.size(); k++)
            {
                if( (L[i].count(k) > 0) && (L[j].count(k) > 0) )
                    if( m->get(X[i], X[k]) < (m->get(X[j], X[k]) - 16.0*epsilon) )
                        return true;
            }
        }

        return false;
    }

    bool twinClusters( std::vector<unsigned int > &clusterOfPsy, unsigned int a, unsigned int b )
    {
        unsigned int ca = clusterOfPsy[a];
        unsigned int cb = clusterOfPsy[b];
        if(ca < cb)
            return ( (ca+1) == cb);
        return ( (cb+1) == ca);
    }

    void calculateDx( std::vector<double> &d, std::vector<unsigned int> &X, Matrix *m )
    {
        d.resize(X.size());
        for (unsigned int i = 0; i < X.size(); i++)
        {
            double min, max;

            if( i == 0 )
                min = max = m->get(X[i] ,X[X.size()-1]);
            else
                min = max = m->get(X[i] ,X[0]);

            for (unsigned int j = 0; j < X.size(); j++)
            {
                if( j != i )
                {
                    if( m->get(X[i], X[j]) < min )
                        min = m->get(X[i], X[j]);
                    if( m->get(X[i], X[j]) > max )
                        max = m->get(X[i], X[j]);
                }
            }
            d[i] = (min + max) / 2.0;
        }
    }
}

